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We investigate the rheology of granular materials via molecular dynamics simulations of homoge- 
neous, simple shear flows of soft, frictional, noncohesive spheres. In agreement with previous results 
for frictionless particles, we observe three flow regimes existing in different domains of particle vol- 
ume fraction and shear rate, with all stress data collapsing upon scaling by powers of the distance to 
the jamming point. Though this jamming point is a function of the interparticle friction coefficient, 
the relation between pressure and strain rate at this point is found to be independent of friction. 
We also propose a rheological model that blends the asymptotic relations in each regime to obtain 
a general description for these flows. Finally, we show that departure from inertial number scalings 
is a direct result of particle softness, with a dimensionless shear rate characterizing the transition. 

PACS numbers: 45.70.-n, 47.57.Gc, 64.60.F-, 64.70. ps, 83.10.Gr, 83.80.Fg 



I. INTRODUCTION 



spheres i and j are 



Flows of granular matter occur in numerous geophysi- 
cal and industrial processes and, as such, have garnered 
the attention of researchers for many years. Early ef- 
forts to describe these flows focused on either dilute flows 
(where kinetic theories [IH1] apply and which belong to 
the inertial regime) or very dense, slow flows (or quasi- 
static flows, for which plasticity models E] can be 
used). However, attention has turned recently to the in- 
terface between these two regimes in the context of a jam- 
ming transition, proposed to occur in granular and other 
soft matter [7] . Of particular interest are several works 
that find a critical rheology around this transition in flows 
of frictionless, soft spheres [BHTT] and disks [TU [T3j ■ Fur- 
thermore, they find scalings for the mean normal and 
shear stresses with respect to volume fraction that apply 
over a wide range of volume fractions and shear rates. 
Granular materials, though, are typically considered stiff, 
frictional materials, and to date there has been little work 
on identifying a critical rheology [Ml [15] for such matter 
despite significant progress in understanding their static 
jamming behavior 16-19 . In this paper, we investigate 
the rheology of frictional granular matter about the jam- 
ming transition and discuss the construction of a rheo- 
logical model for flows in the quasi-static, inertial, and 
intermediate (i.e. critical) regimes. 
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for overlap distance dij, particle diameter d, spring stiff- 
ness constants k n and k t , viscous damping constants 7„ 



and 7t, effective mass m B s = mirrij / (rrii 



for par- 



ticle masses mi and rrij, relative particle velocity com- 
ponents v™- and vjj, and elastic shear displacement u' -. 
A linear spring-dashpot (LSD) model is chosen by set- 
ting the function f(x) = 1, while a Hertzian model 
is set by f(x) = ^/x; the LSD model will be used 
throughout this paper except where noted explicitly. By 
Newton's Third Law, particle j experiences the force 
Fji — — Fy. Particle sliding occurs when the Coulomb 
criterion |F*-| < /i|F™| is not satisfied for particle fric- 
tion coefficient /i. Additionally, after setting k t /k n = 2/7 
and 7t = 0, we set 7„ such that the restitution coefficient 

e =exp^— 7„7r/ \J 4k n /m c fi — 7 2 ^ = 0.7 in the LSD case. 

Using the above contact model, assemblies of about 
2000 particles in a periodic box are subjected to homo- 
geneous steady simple shear at a shear rate 7 via the 
Lees-Edwards boundary condition |22j . The box size, 
and hence the solids volume fraction <f>, are kept con- 
stant for each simulation. The macroscopic stress tensor 
is calculated as 



II. SIMULATION METHODS 

We perform computer simulations using a package of 
the discrete element method (DEM) [20] implemented 
in the molecular dynamics package LAMMPS [3T]. In 
DEM, particles interact only via repulsive, finite-range 
contact forces. We employ a spring-dashpot model, for 
which the normal and tangential forces on a spherical 
particle i resulting from the contact of two identical 
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(3) 



where V is the box volume, ry is the center-to-center 
contact vector from particle j to particle i, and is 
the particle velocity relative to its mean streaming ve- 
locity; from this result, an ensemble- averaged pressure 
V = {°~xx + o~yy + 0zz)/3 and shear stress r = a xz can be 
extracted. All macroscopic quantities will be presented 
in dimensionless form, scaled by some combination of 
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TABLE I. Estimates of the critical volume fraction cj> c for dif- 
ferent cases of the interparticle friction coefficient. The value 
of 4> c for the frictionless case agrees with the experimentally 
determined result of Nordstrom [10] . 
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the particle diameter d, stiffness k — k n , and solid den- 
sity p s . Since particles are assumed to overlap without 
deformation, we ensure that the average overlap is small 
(i.e. S/dxpd/k < 0.07). 

III. FLOW REGIMES 

We performed a series of simple shear simulations over 
a range of shear rates and volume fractions reaching into 
all three flow regimes and for several particle friction coef- 
ficients between and 1. Figure [I] shows the scaled pres- 
sure pd/k versus the scaled shear rate j — jd/ \Jkj (p s d) 
at various volume fractions for (a) p — 0.5 and (b) 
p = 0.1. At low shear rates, there is an observed separa- 
trix occurring at a critical volume fraction </> c , which we 
identify as the jamming point; stresses scale quadratically 
with shear rate below <p c but show no rate dependence 
above it. These two bands correspond to the inertial 
and quasi-static regimes, respectively. As shear rate in- 
creases, the quasi-static and inertial isochores approach 
a shared asymptote characteristic of the critical point in 
which dependence on the volume fraction vanishes; this 
region corresponds to the intermediate regime. Interest- 
ingly, the intermediate asymptote appears to be indepen- 
dent of the friction coefficient, in contrast to results at 
lower shear rates and despite the fact that <f> c = 4> c {p)- 
Values of <fi c for different cases of p are presented in Ta- 
blelU It should be noted that these critical values inferred 
from dynamical behavior of sheared systems are unique 
for each case of p and hence may differ from the jam- 
ming points of static packings, which are not unique and 
depend on the compactivity |17j . 

A better understanding of the regime transitions can 
be gained by constructing a regime map, or "phase dia- 
gram," from the slopes of the curves in Figure [T] Such 
a map is shown in Figure [2j The intermediate regime is 
observed to lie in a window centered around 4> = 4>c, and 
the width of this window is dependent on the value of the 
dimensionless shear rate. This feature has important im- 
plications for the modeling of dense granular flows. The 
large stiffness of granular materials such as sand or glass 
beads has been used to justify the modeling of granu- 
lar particles as (infinitely) hard spheres. For such par- 
ticles, dimensional analysis requires the traditional Bag- 
nold scaling of the stresses (i.e. p, t ~ 7 2 ), thereby ren- 
dering the intermediate and quasi-static regimes impos- 
sible. This picture is consistent with the vanishing of 
the intermediate regime observed in Figure [2] as k — > oo. 



However, real granular materials do nevertheless have a 
finite stiffness. Therefore, in the context of building a 
general rheological model for granular flows, it is prefer- 
able to choose a framework that include all three regimes. 

Another important observation from Figure [2] is the 
smoothness of the transitions between the regimes. This 
feature suggests that purely quasi-static, inertial, or in- 
termediate flow is achieved only in certain limits. As 
7 — > 0, we see quasi-static flow for (f> > <j> c , inertial flow 
for (j> < cf> c , and intermediate flow at <f> = <p c . We also 
see intermediate flow as 7 — > 00 for all volume fractions 
over the wide range examined in this study. The smooth 
transitions also suggest that the rheology at a particu- 
lar (7, 4>) is a composite of contributions from low-7 and 
high-7 behaviors, and this notion will play a large role in 
our construction of a rheological model. 



IV. CRITICAL VOLUME FRACTION cf> c 

Because <p c plays such an important role in governing 
the rheology in each of the three flow regimes, accurate 
estimation of its value for each case of p is required for the 
construction of a valid rheological model. However, this 
task is made difficult by fluctuations of our measurements 
in time t. We observe a propensity for assemblies near <j> c 
to form and break force chains intermittently during the 
shearing process, resulting in stress fluctuations of several 
orders of magnitude as seen in Figure (3}i. Though fluc- 
tuations occur at all volume fractions, their size relative 
to the mean is markedly large near the critical point. In 
Figure^ the standard deviation a p = \J (p 2 {t)) — (p(t)) 2 
of the pressure, when scaled by the time-averaged pres- 
sure p, exhibits a spike centered slightly under <fi c . This 
phenomenon increases the potential error in the time- 
averaged stress values near the critical point, thereby 
limiting the precision of our <f> c estimates to within about 
±0.001. 

Additionally, though <j) c is certainly an important 
quantity, it is not necessarily the only or even the most 
influential parameter describing the jamming transition. 
The fact that stress can vary significantly at a constant 
volume fraction indicates that, while is a useful predic- 
tor of time-averaged stresses, other state variables may 
be more suitable for predicting instantaneous stresses. 
This quality has been observed previously with the co- 
ordination number Z(t), for example, which was shown 
to exhibit a one-to-one correspondence with p(t) in the 
quasi-static regime [23J. Indeed, we observe that p(t) and 
Z(t) exhibit the same qualitative evolution in time (Fig- 
ure [3]:) and a similar 0-dependence in their fluctuations 
(Figure [3ji). Here we define Z(t) = 2N c (t)/N for N c 
contacts occurring in the iV-particle assembly. The p(t)- 
Z(t) relationship suggests that the critical point is better 
defined by some critical coordination number Z c . How- 
ever, because our goal is the construction of a steady-state 
rheological model, it is convenient to ignore all dynamics 
and assume that Z c and <j) c correspond to the same con- 
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FIG. 1. Dimensionless pressure vs. dimensionless shear rate for various volume fractions with (a) n = 0.5 and (b) fi — 0.1. In 
both cases, three flow regimes are observed, each with the scalings p ~ -y m : a quasi-static regime with m = 0, an inertial regime 
with m = 2, and an intermediate regime with m ~ 1/2. At low 7, a critical volume fraction (f> c separates the quasi-static and 
inertial regimes; values of </> c (a0 are given in Table [I] 
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FIG. 2. Regime map for p = 0.5, with volume fraction vs. di- 
mensionless shear rate. The intermediate regime occurs only 
at cf> c in the limit 7 — > but emerges from this point to en- 
compass all volume fractions as 7 — > 00. 



ditions. We therefore proceed with <j> c as the definition 
of the critical point for our model. 

In addition to being a function of /it, the critical point 
has also been proposed to change with the restitution 
coefficient e [23], and such a <fr c (e) has been used in a 
kinetic theory for frictional particles [25 . However, our 
DEM results do not support this conclusion, especially 
for frictional particles. As seen in Figure[4^-b for /1 = 0.5 
and p — 0.1, the spike in the pressure fluctuations oc- 
curs at the same volume fraction for a given ll regard- 
less of the value of e, suggesting that <f) c = 4>c(p) only. 
Even for the frictionless case (Figure |4j;) , where fluctua- 
tions tend occur over a wider range of volume fractions, 
there is no clear trend in the peak towards lower <j>. One 
possible reason for the discrepancy is the methods used 
for determining cj> c . Because hard-sphere codes, used in 
Ref. [23], prohibit particle overlaps, they are unable to 
simulate sheared particle systems near or above <j> c |26j . 



This shortcoming limits the performable simulations to 
one side of (f) c , thus requiring the critical point to be 
estimated via extrapolation. Furthermore, while hard- 
sphere methods treat collisions as binary interactions, 
entrance into the quasi-static regime coincides with the 
development of multi-body interactions that persist even 
in the hard-sphere limit [27] . This conflict may render 
even-driven algorithms less accurate at resolving colli- 
sions upon approaching 4> c and perhaps result in an erro- 
neous estimation of the value of <f> c . Soft-sphere DEM, on 
the other hand, enables us to resolve multi-body contacts 
and simulate shear flows at any volume fraction on either 
side of <j) c , thereby allowing us to interpolate the value of 
4> c - For these reasons, we expect the latter approach to 
provide more accurate 4> c estimates. 



V. PRESSURE SCALINGS AND REGIME 
BLENDING 

It has been demonstrated in experimental [10] and 
computational [5] HI H2] studies of frictionless particles 
that stress data will collapse onto two curves (one above 
4> c and one below) upon scaling the stresses and shear 
rate by powers of \<t>— <f» c \, the distance to jamming. This 
idea is consistent with several models of the radial distri- 
bution function, used in kinetic theories for the inertial 
regime, that diverge at close packing [2"5H3"U] . Such a col- 
lapse can be achieved for frictional particles as well, as 
shown in Figure [5] with 

p*=p/\cf>-cj> c \ a %=^/\4>-^ c \ b (4) 

and constitutive exponents a and b. This result for fric- 
tional disks was also found independently in Ref. [15] . 
From the collapse it is clear that an asymptotic power- 
law relationship between stress and shear rate exists for 
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FIG. 3. Characteristics of the critical point for fi = 0.5 and 7 = 3.2 x 10~ 5 . (a) Pressure fluctuations in time are observed 
to become larger near the critical point, (b) The standard deviation of pressure, scaled by the mean pressure, exhibits a 
spike at <f> c . (c-d) The coordination number fluctuations are similar to those of the pressure in terms of both dynamics and 
^-dependence. 
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FIG. 4. Effect of changing the restitution coefficient on the pressure fluctuations for 7 = 3.2 x 10~ 5 . Dotted lines demarcate 
the critical point cf> c . For (a) /n = 0.5 and (b) fi — 0.1, the location of the spike in pressure fluctuations is independent of e. (c) 
Even for the frictionless case, there is little evidence to suggest 4> c = 4> c {e). 
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each flow regime j, and we can write the form of each 
asymptote as 



TABLE II. Values of model constants 



Pj 



10 ~ 0c 



7 



10- 4>c 



(5) 



where ttiqs = 0, mi ncrt = 2, and mint = m * ■ The ex- 
ponents a and b can be fitted from the DEM data, but 
the values are sensitive to the choice of C used [9] and 
hence should be chosen with care. Our inertial regime 
data suggest that pinert ~ \4> ~ 0c |~ 2 , which is consis- 
tent with previous results [31] : quasi-static regime data 
reveal that pqs ~ \4> — 0c| 2 ^ 3 ; and, as noted earlier, 
Pint ~ 10 — 0d°- These trends lead us to set a = 2/3, 
b = 4/3, and m* = 1/2. The m* value is consistent 
with our fits of the intermediate asymptote (Table [I]) and 
with experimental results [TO] [IT], and it is similar to 
other values proposed for frictionless particles using the 
linear spring-dashpot model [5J H] . The value of a used 
in Ref. [23] (a = 1), though different, still yields a decent 
collapse. However, in that work, (j) c is determined by ex- 
trapolation from the quasi-static regime, while here we 
interpolate it from quasi-static and inertial regime data 
and furthermore verify it with stress fluctuation data, as 
described in Section |IV| hence we believe our current C 
values and the resulting a value to be more accurate. We 
also point out that the above scaling exponents depend 
on the contact model used [8j [32] . Based on a small set of 
simple shear simulations with a Hertzian contact model, 
we observe the values of a f=s 1 and m* w 3/4 to be 
larger than in the LSD case by a factor of 3/2, which is 
consistent with previous results for static, jammed sys- 
tems [16] . The value of b, however, remains the same for 
both contact models; note that in both cases a = bra* in 
order to satisfy the functional forms implied by the col- 
lapse. The resulting collapse for the Hertzian particles is 
shown in Figure [6j 

Though the individual regime limits can be described 
using Eq. [5] the transitions between them have yet to be 
modeled. To this end, we employ a blending function B 
of the form 



B(y 1 ,y2) = (y? + y%) 1/ * 



(6) 



with w > yielding an additive blend for the quasi-static- 
to-intermediate transition and w < providing a har- 
monic blend for the inertial-to-intermediate transition. 
Figure [5] demonstrates the use of Eq. [6] with the asymp- 
totic forms of Eq. [5] and w = ±1. The blended model 
is able to capture the pressure behavior continuously in 
shear rate for all three regime limits as well as the tran- 
sitions; moreover, it does so without defining the stresses 
in piecewise fashion over arbitrary shear rate domains. 
Notably, it also predicts the narrowing intermediate win- 
dow around = C in the limit of zero shear rate, as 
the quasi-static and inertial contributions to the stress 
become small near the jamming point. The general form 
of the pressure model based on the Hookean-case results 
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can hence be written as 



PQS + Pint 
(Plnert 



for 
for 



> 
< 



(7) 



with the individual regime contributions defined as 

PQsdA = a QS |0-0 c | 2/3 (8) 
Pintd/k = aint7 1/2 (9) 



Plncrtd/k = 



aincrt7 
10- C |' 



(10) 



The pressure at 4> = <f> c can be calculated using either 
blend, since Eqs. [8] and [lO] yield pqs(0 = 4>c) = and 
Pincrt(0 = 4>c) = oo, which both yield p = p int upon 
substitution into Eq. [7j this case is included with the 
quasi-static blending solely for the sake of simplicity. The 
constitutive parameter «qs is a function of fi, while ai nC rt 
and aint are fairly /i-independent. These and other model 
constants are given in Table [TTJ 

There are a few features of Eqs. [7[ [10] that are worth 
noting. Firstly, for systems above the critical volume 
fraction, the blending function yields a model of Herschel- 
Bulkley form, which has been shown previously to cap- 
ture the shear stress of soft-sphere systems [TOl QT] . Ad- 
ditionally, the individual regime contributions are con- 
sistent with some known scalings. For example, the 
quasi-static pressure is proportional to the particle stiff- 
ness [H|, while pincrt = ai nelt P s{i d) 2 /\cj) - C | 2 rightly 
exhibits no dependence on k [TJ-[4] . Finally, the viability 
of the 0-scaling in Eq.[l0]for all fi values suggests that the 
0c = 0c (m) formulation could be a simple but effective 
step in improving current kinetic theory models. 



VI. DIMENSIONLESS GROUPS AND STRESS 
RATIO MODEL 

It is possible to construct an analogous model for the 
shear stress as for the pressure, as previous works have 
shown t to exhibit similar scalings with respect to the 
distance to jamming [9j [13] . However, because r and p 
both vary over several orders of magnitude, fitting them 
directly can result in poor predictions of their ratio, i.e. 
the shear stress ratio rj = r/p, which varies over a much 
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FIG. 5. Collapse of pressure vs. shear rate curves from Figure [I] for (a) fj, = 0.5 and (b) fx = 0.1. In both cases, the pressure is 
scaled as p* — p/\<j> — 4> c \ 2 ^ 3 and the shear rate as 7* = "//[</> — (f> c \ 4 ^ 3 . A simple blending function (solid lines) captures regime 
asymptotes and transitions. 




FIG. 6. Collapse of pressure vs. shear rate curves for Hertzian 
particles with [i — 0.5. The volume fractions (and legend) 
are the same as from Figure [ija). Here, p* — p/\(j> — (fee] 1 , 
7* = ~y/\<j> — 4>c 1 4 ^ 3 , and m* « 3/4. Regime asymptotes and 
transitions are captured by the same blending function (solid 
lines) as in the Hookean case. 



narrower range. For this reason, we choose to construct a 
model for r\ and then express the shear stress as r = r)p. 

Some recent, successful rheological models for dense 
granular flows employ a dimensionless parameter called 
the inertial number as the basis for achieving stress col- 
lapses over a range of volume fractions and shear rates 
[3T| 133] 134] . This inertial number / = yd/ \/p/p s is a 
ratio of the timescales of shear deformation and particle 
rearrangement, and the physics of granular flows of hard 
particles is said to be determined by the competition of 
these two mechanisms. When the particles have a finite 
stiffness, however, the binary collision time is nonzero 
and therefore presents yet another important timescale. 
With this point in mind, we note that the dimension- 
less shear rate 7 identified earlier is in fact the ratio of 
the binary collision time to the macroscopic deformation 
time |14j , and we show here that it can be used along with 



the inertial number to characterize soft particle rheology. 

In Figure [7] we plot the stress ratio versus the inertial 
number for fi — 0.5. For the densest systems, i.e. fol- 
low /, r\ exhibits a constant-value asymptote that we 
identify as the yield stress ratio r/ s = i]s(^)i values of 
rj s for different cases of /x are presented in Table [IT] As 
/ increases, 77 then also increases. These same obser- 
vations were made in previous studies of particles in the 
infinitely-hard limit |3T] G22 El] ■ However, unlike in these 
works, we also observe significant scatter as / becomes 
larger, which we will now show to be a consequence of 
the particle softness. 

Because the inertial number models are designed for 
hard particles, we first limit our analysis to cases in which 
particle softness has little effect, i.e. for small 7. Indeed, 
quasi-static and inertial regime data of rj versus / from 
our DEM simulations collapse onto a single curve, with 
the quasi-static regime occurring for / < 10~ 2 and in- 
ertial regimes occurring for / > 1CP 2 . This collapse is 
seen in the inset of Figure [7^ for fi — 0.5. We model this 
curve as 

fftardOO = V*W + (J /J)/3i + I ' ( U ) 

where 1$, ai, and /?i are parameters dictating the transi- 
tion from quasi-static to inertial flow. This form is sim- 
ilar to that of Jop et al. [34] . Interestingly, the increase 
of 77 from rj s is nearly identical for all cases of /i > 0.1 
(Fig. [7)3). Since the interparticle friction coefficient for 
most real granular materials falls in this range, we con- 
veniently take one set of constitutive parameter values 
as suitable averages for our model; these values are pre- 
sented in Table ITTl 

The form of 77hard presented in Eq. [TT] is not the only 
viable option. Another possibility is a simple power law, 
which can be written as ?7hard(-0 = Vs + ct^I^ 1 . This 
form has been used previously by da Cruz and cowork- 
ers |31j with P[ = 1. A comparison between this form, 
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with f}[ — 1 and a\ = 0.6, and the one in Eq. 11 are 
shown in Figure [7]b. The two models agree closely for all 
values of I < 0.3, with a departure occurring for larger 
I. However, with the inertial number models, we need 
to be concerned only with volume fractions greater than 
the freezing transition 4>f = 0.49 [28], where traditional 
kinetic theories fail [25] 132] • At <fif, the kinetic theory 
of Garzo and Dufty pQ predicts / = 0.83, which is con- 
sistent with our DEM results and beyond which we can 
ignore disparities in the 7?hard predictions between the 
two models. Hence, though we continue with Eq. we 
view both for ms a s being acceptable. 



Though Eq. 



11 



captures low- 7 behavior well, inclusion 
of higher-7 cases reveals a noticeable departure from the 
^hardl-O curve, as seen in Figure [7j Specifically, for a 
given value of /, the value of r\ from an intermediate- 
regime flow is consistently lower than that given by the 
Eq. [TTJ This deviation is a consequence of particle soft- 
ness and, in the context of our regime blending, grows 
in magnitude with the intermediate-regime contribution 
to the pressure. Figure [8^, shows the connection between 
the magnitude of this departure 77 so f t = ??hard — V and 7. 
This softness effect, similarly to 7/hard, can be modeled as 



a 2 



(7„/ 7 )ft> + l 



(12) 



where 70 = 0.1, 0L2 — 0.2, and fa — 1 are constants 
describing the transition to intermediate flow. Finally, 
we can write 



77(/,7) = 7/hardU) - %oft(7)> 



(13) 



and, by plotting 77* = i] — r/ s + 77 so f t vs. / as in Figure [8]d, 
we arrive at a collapse of the stress ratio data from all 
three regimes. 



VII. GENERALIZED CONTINUUM MODEL 



Our rheological model therefore consists of Eqs. [7] 



10 for the pressure and Eqs. 11-13 for the shear stress 



ratio. Though the collapses can generally be improved by 
allowing the fitting parameters to be functions of /i rather 
than constants, the fits are nevertheless fairly good and 
hence justify the use of simpler forms. 



While this model was developed for simple shear flows, 
it can be recast to handle general deformation types as 
done in Ref. [23] ■ First, we note that the strain rate 
tensor for simple shear flows is D = -^j(e x e z + e z e x ) 
where e.; are the unit vectors in the i direction. This 
expression can be rearranged to yield 

7 = 2|D|, (14) 



where |D| = y |D T : D is the modulus of D, and D is 
taken to correspond to general deformation types. Fi- 
nally, we write the stress tensor as 

a = p(I-r,S) (15) 

where p and 77 are given by our model, I is the identity 
tensor, and S = S/|D| with S = D |tr(D). Eqs. : 
and [15] allow our rheological model to handle flows m 
more complex geometries as are commonly found in real 
flow scenarios. 



VIII. SUMMARY 

We have investigated shear flows of dense frictional 
granular materials in all three flow regimes in order to 
gain a better understanding of the scalings within each 
regime and the transitions between them. We find scaling 
relations for the pressure with respect to both shear rate 
and the distance to the jamming point and, for the in- 
termediate regime, observe identical power-law behavior 
for particles with different friction coefficients. Further- 
more, we propose a simple blending function for patch- 
ing each regime's asymptotic form in order to predict 
pressure in between regimes. Finally, we decompose the 
shear stress ratio into contributions from two dimension- 
less shear rates, enabling us to quantify the effect of par- 
ticle softness. These findings establish a framework for a 
global model for steady-state simple shear flows of dense 
granular matter. 
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